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. Abstract 
a ; 

We present a numerical analysis of the validity of classical and generalized hy- 
Q^- drodynamics for Lattice Boltzmann Equation (LBE) and Lattice BGK methods 

Q I in two and three dimensions, as a function of the collision parameters of these 

H i models. Our analysis is based on the wave-number dependence of the evolution 

"-^ . operator. Good ranges of validity are found for BGK models as long as the relax- 

;_j ■ ation time is chosen smaller than or equal to unity. The additional freedom in the 

■ - - ' choice of collision parameters for LBE models does not seem to give significant 

improvement. 



^Permanent address: 

Department of Physics and Centre for the Physics of Materials, McGiU University, Rutherford 
Building, 3600 University Street, Montreal, Quebec H3A 2T8, Canada 



Recently, Lattice-gas Automata (LGA) methods have been developed as a new 
computational approach to fluid dynamics Using purely local Boolean oper- 
ations to represent particle collisions, they have proved to be extremely efficient, 
although, due to the fluctuations inherent in the method, statistical averaging 
is necessary in order to extract information. The Lattice Boltzmann Equation 
(LBE) 1^ and Lattice BGK method |Hl @], by using continuous distribution 
functions, eliminate this statistical noise and offer signiflcant computational ad- 
vantages. 

It is important to deflne ways in which these methods can be optimized. In 
the past, since one of the main applications has been the study of flows at high 
Reynolds number, it has been customary to minimize the viscosity by tuning the 
parameters of the collision operators. However, since applications to low Reynolds 
number flows are becoming increasingly important |Q , other criteria might be 
used. In particular, it is of interest to deflne the spatial scale for which the models 
reproduce hydrodynamics, and how this scale depends upon the parameters of 
the simulation. 

This analysis was flrst considered by Luo et al. 0, and developed by Grosflls 
et al. 1^ and Das and Ernst |0 for the study of LGA's. They show that some of 
the simplest LGA models reproduce classical or even generalised hydrodynamics 
only over very large spatial scales, and point out that observations in the 
literature of "negative viscosities" |TU| can be traced to those scale effects. Since, 
however, there have been no similar investigations of LBE or BGK methods, we 
present in this letter a numerical analysis of the validity of hydrodynamics for 
these methods in two and three dimensions, the latter being an extension of the 
analysis of four- dimensional models based on the FCHC lattice [0. 

In the LBE with enhanced collisions [0, the collision operator is linearized 
around the equilibrium distribution function to give the kinetic equation: 

/,(X + C„, t + 1) = /,(X, t) + j2 [/Mx, t) - rr\^, t)] (1) 

(3=1 

where fla/3 is the linearized collision operator, /a(x, t) is the occupation number of 
velocity at node x and time t, /^"^"''(x, t) is the chosen equilibrium distribution 
function and b the number of velocity directions. By symmetry, the matrix 
element Qa/3 depends only on the angle 9 between the directions a and [3, and 
by convention is given the value ae |TT[]. In two dimensions, for a single speed 
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model based on the hexagonal lattice [Q] (6 velocity directions), there are four 
possible matrix elements, qq, qqq, a^o and aiso while for the four dimensional 
FCHC lattice (24 velocity directions), there are five possible matrix elements, oq, 
ago, clqq, ai2o and aigo- The conditions of mass and momentum conservation then 
reduce the number of independent matrix elements, leaving two and three such 
elements for the two and four dimensional situations respectively. 

Linear stability analysis requires that the eigenvalues of Q be contained in the 
interval (—2,0). This places further restrictions on the matrix elements, which 
can be made explicit by writing the non-zero eigenvalues as [jlll] : 

A = 6(ao + aeo) 
(T = -6(ao + 2a6o) 

in two dimensions, with multiplicities 2 and 1 respectively, and 

A = ao — 2agQ + aigo 
a = 3(ao - ai8o)/2 
7 = 3(ao + 6ago + ai8o)/2 

in four dimensions, with multiplicities 9, 8 and 2 respectively. It is also useful to 
note that the eigenvalue A is linked to the kinematic viscosity of the fluid by [jlT 

d + 2\X 27 

where d is the dimensionality of the simulation. The eigenvalues a in two di- 
mensions and a and 7 in four dimensions control the decay of the so-called ghost 
fields 0, and therefore a preferred choice for their values is —1, forcing a rapid 
decay of these unphysical fields. 

The BGK model (so-called in analogy with the BGK treatment of the Boltz- 
mann equation |T^) is a further simplification of the LEE model, whereby the 
collision part of the kinetic equation is parameterised by a single relaxation pa- 
rameter T such that [|]: 

/,(x + c„t + l) = /,(x,t)-r-i ■ [/,(x,t)-/^'^»'(x,t)]. (2) 

The relaxation parameter r (which linear stability requires to be larger than 1/2, 
1/2 < r < 1 being called subrelaxation and r > 1 over- relaxation) is linked to 
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the kinematic viscosity by the relation 

_ c\2t-1) 
^~ 2{d + 2) ' 

The value r = 1 plays the same role as a = 7 = —1 in the LBE model since it 
relaxes the distribution function /q,(x, t) to its equilibrium /^^"*'(x, t) in a single 
time step. 

To analyse the hydrodynamic behaviour of these models, we use the properties 
of the evolution operator i?(k) 0. Thus, if 0a (x,/!:) is the deviation of /q,(x, t) 
from the equilibrium distribution, 

0a(x,t) = /a(x,t)-/r"(x,t), 
H(k.) can be defined by Fourier transforming the kinetic equation (1) such that 

m 

|0(k,t + l)>=/7(k)|0(k,t)>. 

Here, < Ca|0(k,t) >= 0a(k,t), H{k) = D(k)/J(0), and H{0) = I + n, where / 
is the unit matrix, and the displacement operator -D(k) is the diagonal matrix 
diag[exp(— ik ■ Ci), exp{—ik ■ C2), ■ ■ ■ , exp{—ik ■ c^)]. The eigenvalues of H(k), 
defined from 

i7(k)|V'A(k) >=e^^«|^A(k) > 

then give information about the transport coefficients corresponding to the colli- 
sion matrix Q. 

In the long-wavelength regime {k 0, where k = |k|), two types of modes 
exist: hydrodynamic modes, related to the conservation laws, with i?e[zA(k)] ~ 
O(fc^), and rapidly decaying kinetic modes, with i?e[2;A(k)] < 0, without any phys- 
ical significance. Transport coefficients are related to the hydrodynamic modes 
IQ. In a model without explicit energy conservation in two (four) dimensions three 
(five) such modes exist, two (two) propagating but damped sound modes (A = ±) 
and one (three) diffusive shear modes (A =-L), with /m[z_|_(k)] = 0. The real part 
i?e[2A(k)] represents damping, and if the imaginary part /m[2;A(k)] = ±Cs(k)k 
is nonvanishing, the mode propagates with speed Cs(k) The wave- vector 
dependent kinematic shear viscosity is defined as ||^ 



3 



while the sound damping constant is defined as 



r(k) = -Re[z±ik)]/k^. 



In classical hydrodynamics (A; — > 0), the transport coefficients are fc-independent 
by definition. However, when this situation does not hold, but the hydrodynamic 
modes are still clearly separated from the kinetic modes, one can speak of a 
generalized hydrodynamic regime ^ , with transport coefficients which are slowly 
varying functions of k. In lattice-based models, the transport coefficients might 
also depend on the direction of the wave vector, k, refiecting anisotropics due to 
the symmetry of the lattice. By computing the transport coefficients through the 
spectral analysis of the evolution operator, and looking at their k-dependence, 
one can judge the range of validity of the classical and the generalized hydro- 
dynamic regime. In a previous analysis 0, for the simplest Lattice-gas FHP-I 
model with a density of p = 1.8, generalized hydrodynamics were shown to be 
valid up to ~ 0.4 for certain directions of k (Figure 3 of reference 0). 

In order to allow for an analysis of H(k) for the Lattice BGK model, we 
require an effective collision matrix. It is easy to verify that, with 



Eq. (1) reduces to Eq. (2), conservation of both mass and momentum is satisfied 
and all the non-zero eigenvalues of Q are equal to — 1/r. This matrix can therefore 
be employed for the analysis. 

Figure 1(a) shows the real part of a typical spectrum obtained for the two- 
dimensional BGK model on a hexagonal lattice with r = 3/4 and k along the 
X direction (parallel to a lattice vector). One can distinguish the hydrodynamic 
(shear and sound) as well as the kinetic modes (i?e[zA(0)] = In |1 — l/r|). Mixing 
of the two kinds of modes happens at A; ~ 3.0. Figure 1(b) displays —Re[zx(}i)]/k'^ 
for the hydrodynamic modes of the same model, and from this figure we conclude 
that classical hydrodynamics is valid up to ~ 1.2 and generalized hydrody- 
namics up to A; ~ 2.1. This conclusion is supported by Figure 1(c) displaying 
Cs(k) = ±Im[zx{i^)]/k. With k along the y-direction, the ranges are 1.5 and 2.3 
respectively. We have studied the same model for values of r ranging from 0.55 
to 1.5, and find that the range of classical and generalized hydrodynamics, given 
by the behaviour of the real and imaginary parts of zx(k), is essentially the same 
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as for the example given above (i.e. for r = 0.75) as long as r < 1. For r > 1, the 
range rapidly decreases so that, for instance, at r = 4/3, classical hydrodynamics 
is only valid up to /c ~ 0.3, and generalized hydrodynamics up to A; ~ 1.3, for k 
along the x-direction. 

Since the two-dimensional LBE model allows for an adjustment of two inde- 
pendent parameters, one might expect a greater scope for tuning the behaviour 
of -2A(k). However, although the quadratic behaviour of —Re[z\(k)] for hydrody- 
namic modes can extend over a greater range, we find that correspondingly, linear 
behaviour of Jm[zA(k)] is found over a smaller range: the overall range of validity 
of generalized hydrodynamics is scarcely improved compared to the Lattice BGK 
model. Our "optimum range" is very comparable with that for multispeed FHP 
Lattice-gas models 0, where the best results were obtained for the 7-bit FHP-III 
model. 

Analysis for the four-dimensional LBE and Lattice BGK models based on the 
FCHC lattice proceeds similarly. Two and three dimensional data can be obtained 
from these four dimensional models by projecting the lattice onto two or three 
dimensions and defining a reduced collision matrix [Q. Our numerical studies 
show that the spectral behaviour of H(k) constructed with the full matrix or the 
reduced matrices is identical, at least for the physically important hydrodynamic 
modes. We therefore present results only for the full four-dimensional FCHC 
lattice. 

For the BGK model, we have again used values of r ranging from 0.55 to 
1.5. Our findings are that classical hydrodynamics is valid up to /c ~ 1.0 and 
generalized hydrodynamics up to /c ~ 2.0, independent of r, as long as r < 
1. As in two dimensions, the ranges rapidly decrease for r > 1, with, at r = 
4/3, classical hydrodynamics being only valid up to A; ~ 0.3, and generalized 
hydrodynamics up to A; ~ 1.3. These ranges, although smaller or comparable 
to those for the two-dimensional models based on the hexagonal lattice, are still 
considerable: they suggests that generalized hydrodynamics is valid down to a 
spatial scale of around three lattice spacings. 

The four dimensional LBE model allows for an adjustment of three indepen- 
dent parameters. As in the two-dimensional case, the range of validity of gen- 
eralized hydrodynamics is hardly changed compared to the Lattice BGK model, 
but the quadratic behaviour of — -Re[zA(k)] for the hydrodynamic modes can be 
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tuned and extended up to /c ~ 1.5 (for example with A = —'^/'^i (7 — ^ — —1 as 
in Fig.2). In general terms, this seems to correspond to the existence of kinetic 
modes which, at k=0, are less clearly separated from the hydrodynamic regime. 

We thus conclude that the computational advantages of the Lattice BGK 
algorithm are complemented by a significant range of validity for classical and 
generalized hydrodynamics. In two dimensions, using a single speed model on 
a hexagonal lattice, the range is as good as that for LGA models. For two- 
and three-dimensional models based on the FCHC lattice, the range is as good 
or better than that for more general LBE methods. The additional parameters 
available in the LBE method seem to give no further advantage. We expect that 
future applications of the Lattice BGK algorithm will exploit this situation. 
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Captions 



FIG.l. (a) Real part of the spectrum for the two-dimensional BGK model with 
T = 3/4 and k || x. The upper of the two hydrodynamic modes {—RelzxCk)] ~ 
as A; — > 0) is the diffusive shear mode (A =-L), while the lower curve corre- 
sponds to the propagating sound modes (A = ±). (b) The corresponding plot of 
—Re[zx{\s.)]/k^ for the hydrodynamic modes. The upper curve corresponds to the 
kinematic viscosity iy(k.), while the lower curve is the sound damping constant 
r(k). (c) The sound velocity Cs(k) = ±Im[z±{'k)]/k for the same model. 

FIG. 2. (a) Viscosity t/(k) (upper continuous curve) and sound damping con- 
stant r(k) (lower continuous curve) for a four- dimensional LBE model, with 
A = —3/2, (T = 7 = — 1 and k || x. The dashed lines correspond to kinetic modes, 
(b) The sound velocity Cs(k) = ±Im[z±{'k)]/k for the same model. 
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